Surface critical behavior of driven diffusive systems with open boundaries 
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Using field theoretic renormalization group methods we study the critical behavior of a driven 
diffusive system near a boundary perpendicular to the driving force. The boundary acts as a 
particle reservoir which is necessary to maintain the critical particle density in the bulk. The scaling 
behavior of correlation and response functions is governed by a new exponent 771 which is related to 
the anomalous scaling dimension of the chemical potential of the boundary. The new exponent and 
a universal amplitude ratio for the density profile are calculated at first order in e = 5 — d. Some of 
our results are checked by computer simulations. 

PACS: 05.40.+j, 05.70.Fh, 64.60.Ak, 66.30.Dn, 72.70.+m 



I. INTRODUCTION 

In order to study the properties of thermodynamic sys- 
tems far from equilibrium physicists have been looking for 
simple models which capture the main features of non- 
equilibrium phenomena. Driven diffusive systems (DDS) 
introduced by Katz et. al. to model fast ionic conduc- 
tors are characterized by a particles conserving dynam- 
ics and a stationary state which does not satisfy detailed 
balance. Their study has led to the discovery of the con- 
nection between the validity of a conservation law and 
the existence of long-range spatial correlations in non- 
equilibrium steady states. 

A simple microscopic realization of DDS is an Ising 
lattice gas with attractive nearest neighbor interaction 
and an external driving force E which prefers particle 
jumps in the direction parallel to E 0,^. The strength 
of the particle attraction can be varied by a temperature- 
like parameter T. Below a critical value Tc{E) particles 
are separated in the stationary state into regions of high 
and low densities, where the interfaces are oriented par- 
allel to the driving force (for E ^ 0). The phase tran- 
sition at Tc{E) is second order. For an infinite driving 
force particle jumps in the direction antiparallel to E are 
suppressed. In this case the phase transition occurs at 
Tc{oo) « 1.41 rc(0), where T'c(O) is the critical temper- 
ature of the two dimensional Ising model ||^. The criti- 
cal behavior of this system has been extensively studied 
by Monte Carlo simulations and renormalization group 
methods (For a recent review see [Q.) 

In most studies of DDS periodic boundary conditions 
in all directions are imposed to avoid surface effects. In 
this case the particles are driven along a ring or torus. 
In more realistic models particles are fed into the system 
at one side and extracted at the other. The asymmetric 
exclusion model is a DDS with hard core repulsion 
(but without nearest neighbor interaction, i.e. T — 00). 
The density profile in a one-dimensional exclusion model 
with a particle source and a sink was investigated nu- 
merically by Krug |^. His results were later confirmed 
and generalized by exact calculations An impor- 



tant result of these works is the 'maximum current prin- 
ciple' which states the following: If the system is placed 
between two particle reservoirs A and B (with the respec- 
tive densities ca and cb with > cb) and the driving 
force points from A to B, then the bulk density takes 
the value Cmax which maximizes the current j under the 
constraint ca > c> Cb, i-e. 



J (Cmax) 



c{j(c) \CA>C> Cb}- 



(1) 



The maximum current density of an Ising lattice gas with 
attractive particle-particle interaction equals its critical 
density 1/2. The low temperature phase of this system 
with open boundaries in two dimensions has been studied 
by Boal et. al. 0. 

In the present paper field theoretic renormalization 
group methods are employed to investigate the effects 
of open boundaries on DDS at the critical point Tc{E). 
We assume that a plane particle source A perpendicular 
to the driving force is located at the left boundary of the 
system (coordinate z = 0) and impose periodic boundary 
conditions in the transverse directions. The effect of the 
particle source is to suppress density fluctuations and to 
maintain a constant density c^i at z = 0. The particles 
are extracted from the system when they reach a sink B 
located at z = i (L 00). 

It is well known that in physical systems with long 
ranged correlations the influence of a surface extends far 
into the bulk. The critical behavior near a boundary is 
governed by universal scaling laws with critical exponents 
that (in general) cannot be expressed in terms of bulk ex- 
ponents (a review is given in Ref. |l^). The applicabil- 
ity of renormalization group methods to investigate both 
static |ri|-|l5|] and dynamic |]l6|-p^ surface universality 
classes is well established. Especially encouraging are the 
results of Ref. where this technique has been used to 
obtain an approximate profile for one-dimensional DDS 
with open boundaries. It turned out that the profile cal- 
culated by renormalization group improved perturbation 
theory (at one-loop order) was in good agreement with 
the exact result of Ref. 0| . 
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In the next section the semi-infinite extension of the 
field theoretic model for DDS at the critical point (in- 
troduced in 1^) is presented. Above the upper critical 
dimension dc — 5 of this model fluctuations around the 
mean field profile can be treated by naive perturbation 
theory. The mean field profile and Gaussian fl uct uations 
for d > 5 are considered in some detail in Sees. Ill and IV 



since the results of this analysis remain qualitatively valid 
for d < 5. In Sec. ^ the renormalization group is used to 
obtain the scaling behavior of Green functions and the 
density profile below five dimensions. Sec. VIl contains 
a discussion. 



II. THE MODEL 

The analysis in the present paper is based on 
the field theoretic model introduced by Janssen and 
Schmittmann |^ to study the critical behavior of a diffu- 
sive system with a single conserved density subjected to 
an external driving force. The model can be written in 
the form of the continuity equation 



dtS + V_Lj_L + dnjn 



0, 



(2) 



where s(r, t) = c(r, i) — c denotes the deviation of the con- 
centration from its average (bulk) value c, and j_L and j|| 
are the respective components of the current perpendicu- 
lar and parallel to the driving force. The explicit expres- 
sion for the current can be motivated by the following 
symmetry requirements [Q: 

(i) Isotropy with respect to the {d — l)-dimensional 
transverse subspace, 

(ii) invariance of the equation of motion under reversal 
of the driving force {E ^ —E) and particle-hole 
exchange ('charge conjugation', s ^ — s), 

(iii) invariance under force reversal and reflection in r|| 
(the coordinate parallel to the force). 

A continuum model satisfying (i)-(iii) describes, for in- 
stance, the long time and large distance behavior of a 
driven Ising lattice gas at its critical density 1/2. Since 
the current in this system is at its maximum value for 
half filling one may use the maximum current principle 
(in a system with open boundaries) to adjust the bulk 
density to the critical value. 

Keeping only terms which are consistent with the 
above conditions (i)-(iii) and relevant or marginal in the 
renormalization group sense the current may be written 
(upon rescaling of s) in the form 



j±=^± [-X{ts- A^s) + C], 
j\\ = E {aQ+ 0-25^) - A/99||S, 



(3a) 
(3b) 



where C is a Gaussian random force with the correlation 



The third order derivative in Eq. (^^ has to be kept 
because the coefficient r of the first order derivative 
vanishes at the critical point (transverse phase transi- 
tion Q). The coefficient cto niay be interpreted as the 
conductivity of the system at the maximum current den- 
sity c. Deviations from c due to fluctuations or a non- 
constant density proflle decrease the current. This effect 
is modeled by the term Ea2S^ in Eq. (^). The coef- 
ficient r measures the deviation of the temperature pa- 
rameter from its critical value and p takes into account 
the anisotropy of the diffusion constant. Even if the diffu- 
sion constant is isotropic (p = 1) in the original Langevin 
equation it becomes anisotropic under coarse graining. 

For the subsequent field theoretic analysis it is con- 
venient to recast the model in the form of the dynamic 
functional M^M' 



Jb[s, s] 



J dt Jyd'^r\^Sdts + X[{A^S){A^s) 

+ r(ViS)(Vi.s)+p(a||.5)(5||s) 
+ ^g{dfs)s^ - hdfs - {V ^s)^]Y 



(4) 



where Xg ^ —Ea2, and s is a Martin-Siggia-Rose re- 
sponse field 1^^. While the functional (^) allows us to 
calculate response and correlation functions for an infi- 
nite bulk system the influence of the boundaries has to 
be modeled by additional surface contributions. If the 
boundary is perpendicular to the driving force the region 
of integration in Eq. (^ is the half space V = {(r_L, z) \ 
rj^ G R'*"-'^, z > 0}, and the surface is defined by z = 0. 
Omitting irrelevant and redundant terms Jl9[ the surface 
functional reads 

Js,[s, s\= J dt Jgy d'^^^r±x(^css — cs^ — C2sA±s 

+ \gs~ss^ - hsS^ (5) 

Response and correlation functions can now be calcu- 
lated by functional integration with the weight exp(— »/), 
where J = Jh + Js- 



III. EQUATION OF MOTION AND MEAN FIELD 
APPROXIMATION 

An exact equation for the stationary profile ^{z) — 
{s{r,t)) follows from the invariance of the generating 
functional 

Z[J, J; Ji, Ji] = y ^[*' exp(-Jb[s, s] - Ji[ss, Ss] 

+ (J, S) + (J, s) + (Ji, S,) + (Ji, s,)) (6) 



under an infinitesimal shift of the field s. In Eq. (|^) 
Ss{r±) — s{r±, 0) denotes the field at the surface and the 
abbreviations 



(C(r,t)C(r',0> -2A<5(r-r')<5(^-0■ 



(3c) 
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(J, S) = I dt I (frJs 



and 



dt I d'^~^rj_JiSs 
dv 



(7) 



(8) 



have been used. The invariance of Z[J , J; Ji, Ji] impUes 
the equation of motion 

dts + A[(A_L - r)A_LS - d\\ {pd\\s + ]^gs^) + 2A_lS] = J 

(9) 

which holds after insertion into averages. The invariance 
of Z[J, J; Ji, Ji] under a shift of the surface field Sg leads 
to the equation of motion 

- pdnS - ^—^sl + CSs - 2css 

-C2Aj_Ss - h, + h=jJi (10) 

which fixes the boundary condition. Taking the average 
on both sides of (^) for vanishing sources J = J = Ji = 
Ji = one obtains 

[pdMz) + Ig^zf + C{z))] = (11) 

or, since 'I'buik = due to the definition of s, 

$'(z) + |- ($(z)2 + C{z) - Cbuik) = 0. (12) 
2p 

The function C{z) = {[s{r±, z; t) — $(z)]^) describes den- 
sity fluctuations at the distance z from the surface and 
C'buik denotes its value for z — > oo. 

In the mean field approximation one neglects the cor- 
relation function C{z) and obtains for the profile 



$mf(2) = *o l + :7-$o2 
2p 



(13) 



Dimensional analysis shows that the momentum dimen- 
sion of the coupling coefficient g is given by [g] = (5 — 
d)/2, and the mean field approximation breaks down be- 
low the upper critical dimension d^ — b. For d > 5 
corrections to the mean field profile can be obtained by 
nai've perturbation theory. At lowest order it is sufficient 
to calculate the perturbation C{z) — Cbuik in Eq. (|l2|) by 
a Gaussian approximation. 



IV. CORRECTIONS TO THE MEAN FIELD 
PROFILE FOR D > 5 

In the simplest case, $„if(2;) = = hg = 0, the 
Fourier transform (with respect to rj^ and t) of the Gaus- 
sian propagator 



G{r±,z,z';t) = (s(ri, z; t)S(0, z'; 0))o 
is given by [ p^ 

1 



(14) 



qi,cj(^; ^ ) 
+ 



-k\z-z\1 ^ 



with 



1/2 



(15) 



(16) 



The parameter c occurring in the surface functional J7i 
and in the propagator describes (for c > 0) the suppres- 
sion of density fluctuations by the particle reservoir at 
the boundary. Since the momentum dimension of c is 
one the asymptotic scaling behavior is governed by the 
fixed point = cx). At this fixed point the fields s and 
s satisfy the Dirichlet boundary conditions Ss — Sg = 0. 

The Fourier transform of the Gaussian correlation 
function 



C{r±,z,z';t) = (s(r_L,z;t)s(O,z';0))o 



(17) 



at the Dirichlet fixed point c^, — oo can be derived from 
the Gaussian part of the dynamic functional. This cal- 
culation yields 



Cqj_,i^iz, z) = 2\q\ 
2Xq 



dz Cqj^^^(z,Z )G'qj^^— a;(^ :^) 



3'[Gq^^i^(z, z')]. 



(18) 



where 3[- • •] denotes the imaginary part. The equal-time 
correlation function at the point (rj^,z) is given by 



C{z) 



>{z,z) 



= Gbuik-^(87rz/Vp)-('^-^)/^ (19) 



with Gbuik ^ ^'''^^/y/P (where A is a cut-off wave num- 
ber). 

We can now use Eqs. (|i^ ) and ( [l9| ) to compute the 
fluctuation correction to the constant mean fleld proflle 
*i'mf(^) = 0. At flrst order in g we get 



<I>W(2 



5(8^) 



-(d-l)/2 



2(rf~3)p 



-{d-3)/2 



(20) 



One can easily check by dimensional analysis that 
higher order corrections to the proflle decay as 
$[2"+il(z)/$[il(z) - 2-"('i-5)/2 (up to cut-off dependent 
terms which may change the amplitude of the leading 
term proportional to z^^"^^^^^^). 

In the limit c, /is ~* oo (with hs/c =: hi fixed) the 
boundary value of the mean field profile is given by $0 = 



3 



hi. This follows from Eqs. (|l^) and (|T|). For $o > 
the mean field profile decays asymptotically as $mf(-2) — 
2p/{gz), and the fluctuation correction ^ 2;-(<^-3)/2 can 
be neglected for z — > cxa (d > 5). However, if hi is positive 
but small the profile $(z) is negative for z < (, where 
C — C^ii) is a crossover length. The dependence of the 
profile on the boundary chemical potential hi is depicted 
qualitatively in Fig. The crossover length ( tends to 
infinity for hi 0"*". To estimate C for small hi we 
equate the mean field profile $mf (C) with the fluctuation 
term ^-('^-3)/2 ^^^^ obtain 



2/{d-3) 



(21) 



for hi ^ 0+, d > 5. 

In the language of semi-inflnite magnetic systems the 
case hi = oo corresponds to the normal transition and 
/ii = is the ordinary point. A length scale similar 
to C has already been found in magnetic systems near 
the ordinary transition |^,^. There the length scale 
C characterizes the magnetization proflle induced by a 
small magnetic surface field. 




FIG. 1. Sketch of the profile ${z) for hi — oo (dotted), 
hi > finite (solid curve), and hi — (dashed). For hi < 
the density in the bulk stays below its critical value indicated 
by the horizontal line. 



V. RENORMALIZATION GROUP ANALYSIS 

A. Renormalization 

The naive perturbation theory described in the pre- 
vious section breaks down below the upper critical di- 
mension c?c = 5. The renormalization group allows us to 
improve the perturbation expansion by a partial resum- 
mation. 

Since the individual terms of the perturbation series 
contain for d = dc ultraviolet divergent integrals a regu- 
larization prescription is required to obtain well-defined 
expressions for the otherwise infinite integrals. Here we 



use the dimensional regularization method (analytic con- 
tinuation of the integrals as functions of d). The re- 
maining poles in e — dc — d are then absorbed into 
reparamctrizations of the coupling coefficients and the 
fields. In the field theory for the bulk model (without 
a surface) a renormalization of the parameter p = Zpn 
is sufficient to cancel the ultraviolet divergences at every 
order of the perturbation theory At one-loop order 
the renormalization factor is given by 



l-- + 0{u'), 



A 2 -3/2 - 
Pr P 



(22) 



where p is an external momentum scale and the geometri- 
cal factor = (3/4)(47r)-''/2F((3-e)/2)F(l-he/2)/F(2- 
e/2) has been introduced for convenience. (The index 'i?' 
indicates renormalized quantities.) 

In order to investigate the scaling behavior of response 
functions near the boundary one has to calculate Green 
functions with insertions of the surface response field 
Ss- Since the Gaussian propagator ( |l5| ) vanishes at the 
Dirichlet fixed point = cxd it is necessary togo to higher 
orders in to obtain non- trivial resufis At 
first order in the propagator becomes (for z' = 0) 



+ . . 



(23) 



where GqJ^uj{z, z') denotes the propagator ( p^ ) for c = 
oo. This shows that the leading order terms in an expan- 
sion in powers of c^^ can be studied in the framework 
of a field theory with Dirichlet boundary conditions after 
replacing in expectation values 



^pdnS. 



(24) 



Analogously insertions of the surface field Sg have to be 
replaced (at leading order) by the the normal derivative 

Since a boundary breaks the translational invariance 
of the system it gives rise to new divergences in the per- 
turbation series which are located at the surface [i.e., 
proportional to S{z))]. These surface divergences have to 
be subtracted by appropriate counter terms added to the 
dynamic functional J^. In the appendix it is shown that 
the required counter terms have the form 

Jhct[S,s]= J dt j^d-'r\[pR{Z -l){dfs){d\\s) 

+p„''^A,gp-'AT^d\^s] (25) 
to remove bulk divergences and 



Jsct [s, s] 



dt / d^-^^Xlp-^^^^A.gp-'KipdlS) 

JdV 

+B{pdns)ss + p-R^ A,gp~'FT{pdns)\ (26) 



to cancel e-poles located at the surface. The renormal- 
ization parameters A, B, F, K are calculated at one-loop 
order with the result 



4 



A = — 



6e 



K 



The first term in jTbct renormalizes the diffusion coeffi- 
cient p. The bulk counter term proportional to d\\ s cor- 
responds to a renormalization of the bulk current p9[ 



h — hi 



-1/2 A 

Pr A^9P 



(28) 



We now show that the surface counter terms propor- 
tional to the operators {pdns)ss and p&f^s can be re- 
placed by a multiplicative renormalization of the surface 
response field p9„s. In order to study single insertions of 
{pdns)ss and pd^s in Feynman diagrams we first connect 
them to the Gaussian response function G^^^uiz, z'), i.e. 
we calculate the Gaussian expectation values 



= g 



FIG. 2. Effect of the vertex Xpd^s in a Feynman diagram. 
The hatched area represents the boundary 2 = 0. Each short 
line perpendicular to a propagator line indicates a derivative 
with respect to z. 



The above analysis shows that the counter terms pro- 
portional to pd^s and {pdns)ss give at one-loop order rise 
to a renormalization of the surface response field 



[pd„s] 



-1/2 



pd„s, 



(33) 



s{z)S{z') / dt' / d''-^x'X{pd„S)ss ) = 0, 



and 



siz){pdnS) / dt' / d''~^x'X{pdnS)s 



1 -../V7^ - {s{z){pdnS)) 



s{z) I d'^-^x'\pdl~s ) = -6{z). 



(29) 



(30) 



(31) 



where 

Z^^''^ = l-B-uK + 0{u^) = 1 - — + 0{u^). (34) 

36 

The relation between the redundant surface couplings 
(here B, K) and Zi can be extended to higher orders 
in M in a similar way as in the case of the ^''^-model . 

The last counter term in J^sct which couples to pdnS 
renormalizes the surface chemical potential 



hi 



Z-"\[hiU-p-^'A,gp-^FT) 



(35) 



In Eqs. ( pO| ) and (^l|) the vertices with normal derivatives 
have to be interpreted as 



\{pdnS)Ss 



lim \{pd\\s{z))ss{z) 

2— >0 



and 



(32a) 



B. Scaling 

With the renormalizations at hand we are in a position 
to determine the scaling behavior of the Green functions 



Xpd, 



lim Xpd'us{z) 



(32b) 



respectively, where the limit z ^ oo has to be taken after 
the averages (• • •)o have been performed pl] , pO[ . 

Equations ( p9| ) and ( ^0| ) show that the counter term 
BX{pdns)ss has an effect only in diagrams in which it is 
connected to the surface response field p9„s. In Green 
functions with an insertion of pdnS it effectively is re- 
places pdnS with (1 — B)pdnS. 

If the argument of s on the l.h.s. of Eq. ( |3l| ) is fixed 
(with z > 0) the average in ( ^l| ) vanishes. However, 
due to the ^-function one obtains a non-zero result if 
Eq. ( |3l] ) is integrated over z. Such an integration oc- 
curs in the calculation of the Feynman diagram shown 
in Fig. ^, where the counter term vertex proportional to 
Xpd^s is connected to the bulk vertex g. Fig. || shows 
that an insertion of the counter term 



Xpjj^^^Aegp "Kpd'is 
in a Feynman diagram has the same effect as the vertex 

uKXpB{dnS)Ss- 



2 



IN N 

Gf;^^({r,x,t}) = iT{~s{v,,t^\{s{v,,t,) 

j=i 

(conn) 

M M \ 

X n PC'„s(xfe,isfe) J|p9„s(x;,is;) \ . (36) 



fc=l 



with M insertions of the surface response field pdnS and 
M insertions of pdnS. In this subsection we omit the 
index 'i?' since only renormalized quantities are used. 

At the critical point t = h ~ the Green functions 
satisfy the renormalization group equation 



p^ + Piu)-^ + au)p-^ 

op ou op 



M 



Gj^^_'j!^^({i", X, t}\ u, p, hi;X, p) 







(37) 



which follows from the independence of the unrenormal- 
ized Green functions from the momentum scale p. The 
Wilson functions are given by 



5 



C(w) = M 



d\np 





= u 



71 (m) = /i 



d In Zi 



h-^c(.)), 



(38a) 
(38b) 
(38c) 



where the derivatives are calculated at fixed bare param- 
eters. 

The renormalization group equation ( ^7| ) can be solved 
by the method of characteristics with the result 



G;^'*^({r, X, t}; u, p, h^; A, p) = X^{1)''/^ 



xGj^;^({r, X, t}; u{l),X{l)p, X,{l)^'^h,-X, fil). (39) 

The functions u{l), X{1), and Xi{l) are solutions of the 
set of ordinary differential equations 

dujl) 
dlnl 
d\nX{l) 
dlnl 
d\nXi{l) _ 
dlnl ~ 

In the limit I — > the scale dependent coupling coefficient 
u{l) tends to the fixed point value u^, = (2/3)e + O(e^) 
and the Green functions display power law scaling with 



= /3(u(/)), 


u(l) 




(40a) 


= C(s(0), 


X{1) 




(40b) 


= 7i(S(0), 




= 1. 


(40c) 



1 



X{1) ~ XJ-^\ r, : 
Xi(O^Xi,J''S m=7i(w.) 



C("*) 



4e 



(41) 



+ 0(e2). (42) 



The amplitudes and ^i,* are not universal. 

Eq. (|3^) can be further simplified by dimensional anal- 
ysis. The momentum dimensions of rj^, r||, t, s, s, and 
hi follow from the form of the functional J' (which has 
to be dimcnsionlcss) and are given by 



s^hi^ p-l/V('^"l)/^ (43) 



respectively. For small I Eq. ( |39D maps the large length 
and time scales of the critical region on scales on which 
Green functions can be calculated perturbatively. Here 
we are especially interested in the profile ^{z) = G^'\{z). 
Choosing for the flow parameter the value 



9 



-1/(2+'?) 







(44) 



we obtain from Eq. (^9|) in conjunction with dimensional 
analysis the scaling form 



with the exponents 



(45) 



d-l + T] 
2(2 + ry) 



11 -d 



and 



1 rf-l+^-m 

9 



2(2 + r;) 



(46) 



(47) 



In (45) a and b are non- universal scale factors while the 
scaling function F is universal. 



C. A universal amplitude ratio 

We know from the discussion of the mean field profile 
and the fluctuation corrections in section IV that $(z, hi) 
is finite and non-zero in both cases hi — oo and hi — 0. 
It therefore makes sense to define the universal amplitude 
ratio 



D — lim 



$(z,0) _ F(0) 



(48) 



A perturbative calculation based on the results of sec- 



tion 



IV yields 
<i>(z,oo) 



.9^(8^)- 



-{d-l)/2 



0(^t^^te). 



f — 



4(d-3)p3/2 V^pJ 

6 



(49) 



Upon application of the renormalization group transfor- 
mation with the choice (^^ for the fiow parameter this 
becomes 



D 



(50) 



In the upper critical dimension d = 5 the solution of 
the fiow equation (40a) reads 



u{l) 



31n(l//) 
This yields for the profile 
$(z,0) 2 



for / 0. 



for z 



oo, 



(51) 



(52) 



$(z,oo) 91n(z/zo) 
where zq is non-universal. 



D. Distant wall corrections 

Until now the profile near a particle source has been in- 
vestigated assuming that the particles are extracted by a 
distant sink located at z — L oo. In computer sim- 
ulations only comparatively small systems can by stud- 
ied and corrections to the profile (^5|) due to the distant 
sink become important. At mean field level the profile 



6 



which satisfies the boundary conditions $(0) 
<i>(L) = —CO is given by 



$mf(z) 



2ttp 



cot 



(t) 



2p 



1 /TTzy 
3 vTJ 



oo and 



(53) 



The powers of {z/L) occurring in this expansion below 
the upper critical dimension can be obtained from a short 
distance expansion (SDE) of the order parameter field 
s{z) for z ^ (3l]-|3|]. The leading term in this SDE 
(with the lowest momentum dimension) is the unit op- 
erator 1. Since Ss — Q due to the Dirichlet boundary 
condition the next-to-leading contribution is the normal 
derivative p9„s. We therefore obtain 



(54) 



The power in front of the normal derivative has been de- 
termined by comparing the anomalous scaling dimensions 
of the individual terms in equation (p4) , 



(55) 



The SDE (b4) implies that the distant particle sink gives 

rise to a correction to the profile proportional to z 2+1 = 
z'^ for z ^ 0, i.e. 



(f>(z) = Aiz' 



l + B 



L 



(56) 



For e = this form is consistent with the mean field 
result (H). 

The amplitudes Ai and B depend on the fixed point 
value of the surface potential, i.e. they take different 
values for /ii = and hi =00. Equation (|5^) shows 
that for hi ^ 00 the (universal) amplitude B is given by 
B = — 7r^/3 + 0(e). In the case hi = Q with the boundary 
conditions $mf(0) = and $(i) = —00 the mean field 
profile reads 



^mliz) 



no /7rz\ 
-—tan — 
gL \2LJ 



2gz 



+ 



(57) 



To determine the amplitude B at leading order in e we 
divide $nif (2) by the semi-infinite profile (Uw and obtain 



^mfjz) ^ 3^. 
$W(z) 2u ^ 



o(„.))(^) 



-e/2 



L 



(58) 



At the fixed point Ui, — (2/3)e 4- O(e^) the amphtude B 
is thus given hy B = dn^/iie) + 0(e°). Note that B is 
of the order 1/e because the semi- infinite profile vanishes 
at zero loop order. 



VI. SIMULATION RESULTS 

In order to check some of the results presented in 
the previous section by computer simulations we use the 
standard Monte Carlo technique with Metropolis spin- 
exchange jump rates on the two-dimensional, driven Ising 
lattice gas with attractive interactions . The driving 
force is effectively infinity, i.e. every attempt of a particle 
to jump in the direction of the driving force is successful 
unless the jump would violate the excluded volume con- 
straint. Jumps in the direction antiparallel to the driving 
force have zero probability. We use the critical value of 
the temperature parameter Tc(oo) = 1.41Tc(0) obtained 
by Leung The particle reservoirs at the boundaries 
are incorporated into the model by a simple change of 
the updating algorithm: Whenever a boundary site is 
involved in an updating step the occupation number of 
this site set equal to a random number X £ {0, 1} which 
is one with propability ca (at the left boundary) or cb 
(at the right boundary). To avoid unwanted correlations 
each realization of X has to be used for only one update. 
In the transverse directions periodic boundary conditions 
are imposed. 




1000 



FIG. 3. Density profile for ca ~ 1.0, cb ~ 0.5, where 
the occupation numbers are represented by the spin variable 
CT = 2n — 1. The statistical error is everywhere smaller than 
±0.006. The solid line is a fit using Eq. @ with Ai = 0.678 
and B = -1.62. 

Fig. H shows the density profile for ca — 1.0 and 
Cb — 0.5. The system size is L|| = 1000 in the direc- 
tion parallel to the driving force and L± = 500. At the 
beginning of each run, an uncorrelated initial state is 
generated where each lattice site is occupied with proba- 
bility 0.5. Then 10^ Monte Carlo steps (per site) are per- 
formed to reach the stationary state. The profile shown 
in Fig. ^ has been obtained by averaging over 2-10^ 
configurations. The amplitudes Ai and B in Eq. ( |5^ ) 
have been determined by a least square fit with the re- 
suh Ai = 0.678 ± 0.004, B = -1.6 ± 0.2. For this fit 
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we have used various subintervals of 3 < z < 50. (The 
statistical error in this range is smaller than 0.002.) For 
larger values of z higher powers in z/L become increas- 
ingly important. We have checked that the above values 
for A\ and B also provide acceptable fits for smaller sys- 
tems = (500,397) and (125,250)]. 



0.2 
0.15 

.1 

0.07 
0.05 

0.03 




500 



Z 



FIG. 4. Double logarithmic plot of the density profile for 
CA = 0.278, 0.280, 0.282, 0.284, 0.286 (from top to bottom) 
and cb ~ 0.5. The spin variable cr = 2n — 1 has been used. 
The broken line corresponds to the power 0.292"^'''^. 




FIG. 5. The density profile for ca = 0.278, 0.279, 0.280, 
0.282, 0.284, 0.286 (from bottom to top) and cg = 0.5. 
The occupation numbers are represented by the spin vari- 
able a = 2n — I. The solid curve is a fit using Eq. (5^ with 
Ai = -0.30 and B = 0.51. 

To determine the amplitude ratio D one first has to 
find the critical value of ca that corresponds to a van- 
ishing surface field, hi =0. Fig. ^ suggests that this 
value is close to ca ~ 0.28. For 0.278 < ca < 0.282 we 
obtained fits consistent with Ai = —0.29 ± 0.02. One 



of these fits is depicted in Fig. |^ together with density 
profiles for various values of ca ■ Each profile is an aver- 
age over 10^ configurations. To determine the amplitude 
B it would be necessary to obtain a more accurate es- 
timate for the critical surface density. The simulation 
result for the amplitude ratio reads D = —0.43 ± 0.03 
which can be compared with our one-loop calculation, 
D « -e/9 = -0.33 for e = 3. 



VII. SUMMARY AND OUTLOOK 

A particle reservoir coupled to the boundary of a driven 
diffusive system maintains the critical density in the bulk 
if the chemical potential of the reservoir is not below a 
critical value. Above this critical value the density pro- 
file (as a function of the distance from the boundary) 
asymptotically approaches the bulk density from above, 
where the decay of the profile follows a power law with an 
exponent cr which can be expressed in terms of the bulk 
exponent rj. At the critical value of the boundary chem- 
ical potential the density tends to its bulk value from 
below. If the chemical potential is close to (but above) 
its critical value the density profile crosses the critical 
density at a macroscopic distance C from the boundary. 
The singular power law dependence of the length scale ^ 
on the boundary chemical potential is characterized by a 
new exponent vi which has been calculated at first order 
in e = 5 — d. 

While in exclusion models without particle-particle at- 
traction the density profile is always a monotonic func- 
tion of the distance from the boundary we have shown 
that in critical DDS stationary profiles can have local 
maximum points. This is due to the density correlations 
in the bulk generated (for d > 1) by the attractive inter- 
action. If the 'temperature' is raised above Tc{E) these 
correlations survive as long as T is finite. Therefore the 
qualitative form of the density profile will not change for 
Tc{E) <T <oo. 

In this paper one out of a multitude of universality 
classes describing various types of DDS has been con- 
sidered. These universality classes differ in the nature 
of the noise (particle conserving or non-conserving) , the 
presence or absence of quenched disorder and the val- 
ues of temperature-like critical parameters . We plan 
to extend the analysis presented here to other universal- 
ity classes. It is straightforward to derive relations simi- 
lar to ( ^ ) between a and the anisotropy exponent r] for 
DDS with quenched disorder. This makes it possible to 
check the field theoretic predictions of Refs. ||3^,^,Q for 
disordered DDS by Monte Carlo simulations of the den- 
sity profile in systems with open boundaries. Note that 
in the presence of quenched disorder periodic boundary 
conditions (in the direction parallel to the driving force) 
may lead to unwanted correlations since the particles are 
subjected to the same randomness after every passage 
through the system. 



8 



In order to obtain a numerical estimate for the surface 
exponent i^i or a more accurate value for the amplitude 
ratio D it is necessary to determine the critical surface 
density more accurately. This is an open problem for 
future simulations. 



ACKNOWLEDGMENTS 

This work has been supported in part by the Son- 
derforschungsbereich 237 [Unordnung und Grofic Fluk- 
tuationen (Disorder and Large Fluctuations)] of the 
Deutsche Forschungsgemeinschaft. 



APPENDIX A: SURFACE DIVERGENCES AT 
ONE-LOOP ORDER 

In order to determine the renormalization constants 
at one-loop order one has to evaluate the ultraviolet di- 
vergent Feynman diagrams shown in Fig. |[ The results 
have to be interpreted in the distribution sense since the 
calculation of Green functions involves integrations over 
the z-coordinates of amputated graphs. 

The Laplace transform of the first diagram in Fig. ^ 
reads 



dze 



x(-^ + | + |Vp.s + 0(.)), (Al) 

where Cq^^^{z, z) is the Gaussian correlator (|lj) at the 
Dirichlet fixed point. The two dimensional Laplace trans- 
form of the second diagram is given by 



(XgY I dze 

10 



dz'e 



{Z, z')dz'Gq^^u,{z, z') 



(A2) 



_xl_ 

2^e 



\s + s' 



0{e) 



Appl ying the inverse Laplace transformation to (Al) 
and (A2) we obtain 



Graph 6(a) = -^A,t"/^(— + —S{z) 
e V^Ao 3 



2 4r, 
+ lVpS'{z) + 0{e)) 



(A3) 



and 



Graph 6(b) = 



■At- 



^/2(2J'(z' I z) 



ls{z')S{z) + 0{e)), (A4) 



where we have introduced the definition 



dz'S\z'\z)f{z') = -f{z). 



(A5) 



(a) 



(b) 



FIG. 6. Ultraviolet divergent Feynman diagrams at 
one-loop order. A line with (without) an arrow represents 
the Gaussian propagator (correlator). The short line perpen- 
dicular to the propagator hne in the diagram (b) indicates a 
derivative with respect to z' . 



The e-poles are cancele d by the counter terms (& 



and (|26|) given in section V A. The values of the coef 



ficients A, B, F, K, and Z at one-loop order follow from 
Eqs. (^) and (^) as 



A = — 
e 

and 



3e 



4 



Z = 1- -. 
e 



K=- (A6) 



(A7) 
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